#include "Itera.h"
#include "Basic_Opera.h"
#include "Result.h"

/**
* @Name: RK4 
  @Brief: 4th-order Runge_Kutta method in t
*/
void RK4()
{
	fftw_complex *C_phi, *C_psi;
	C_phi = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*(M*NUM));    /// rhs
    C_psi = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*(M*NUM));
	Step = 0;

	for(i = 0; i < M*NUM; i++){											  /// start
		C_omega_o[i][0] = C_omega[i][0];
		C_omega_o[i][1] = C_omega[i][1];
		C_rho_o[i][0]   = C_rho[i][0];
		C_rho_o[i][1]   = C_rho[i][1];
	}

	do{																	  /// Write result
		if(Step % 100 == 0){
			Result(Step);
			if(myrank == 0){
				printf("The itera step is %d, writing results...\n", Step);
			}
		}

		for(i = 0; i < M*NUM; i++){  										  /// 1st, get K1
			C_omega[i][0] = C_omega_o[i][0];
			C_omega[i][1] = C_omega_o[i][1];
			C_rho[i][0]   = C_rho_o[i][0];
			C_rho[i][1]   = C_rho_o[i][1];
		}	
		Value_u_v(C_omega, C_u, C_v);										  /// get u, v
		Convolution(C_phi, C_psi);											  /// cal convolution
		for(i = 0; i < M*NUM; i++){
			C_omega_n[i][0] = C_omega_o[i][0] - DH_t*(C_psi[i][0] + C_rho_x[i][0])/6;
			C_omega_n[i][1] = C_omega_o[i][1] - DH_t*(C_psi[i][1] + C_rho_x[i][1])/6;
			C_rho_n[i][0]   = C_rho_o[i][0]   - C_phi[i][0]*DH_t/6;
			C_rho_n[i][1]   = C_rho_o[i][1]   - C_phi[i][1]*DH_t/6;
		}
	
		for(i = 0; i < M*NUM; i++){											  /// 2nd, get K2
			C_omega[i][0]   = C_omega_o[i][0] - DH_t*(C_psi[i][0] + C_rho_x[i][0])/2;
			C_omega[i][1]   = C_omega_o[i][1] - DH_t*(C_psi[i][1] + C_rho_x[i][1])/2;
			C_rho[i][0]     = C_rho_o[i][0]   - C_phi[i][0]*DH_t/2;
			C_rho[i][1]     = C_rho_o[i][1]   - C_phi[i][1]*DH_t/2;
		}
		Value_u_v(C_omega, C_u, C_v);
		Convolution(C_phi, C_psi);
		for(i = 0; i < M*NUM; i++){
			C_omega_n[i][0] = C_omega_n[i][0] - DH_t*(C_psi[i][0] + C_rho_x[i][0])/3;
			C_omega_n[i][1] = C_omega_n[i][1] - DH_t*(C_psi[i][1] + C_rho_x[i][1])/3;
			C_rho_n[i][0]   = C_rho_n[i][0]   - C_phi[i][0]*DH_t/3;
			C_rho_n[i][1]   = C_rho_n[i][1]   - C_phi[i][1]*DH_t/3;
		}

		for(i = 0; i < M*NUM; i++){										       /// 3rd, get K3
			C_omega[i][0]   = C_omega_o[i][0] - DH_t*(C_psi[i][0] + C_rho_x[i][0])/2;
			C_omega[i][1]   = C_omega_o[i][1] - DH_t*(C_psi[i][1] + C_rho_x[i][1])/2;
			C_rho[i][0]     = C_rho_o[i][0]   - C_phi[i][0]*DH_t/2;
			C_rho[i][1]     = C_rho_o[i][1]   - C_phi[i][1]*DH_t/2;
		}
		Value_u_v(C_omega, C_u, C_v);
		Convolution(C_phi, C_psi);
		for(i = 0; i < M*NUM; i++){
			C_omega_n[i][0] = C_omega_n[i][0] - DH_t*(C_psi[i][0] + C_rho_x[i][0])/3;	
			C_omega_n[i][1] = C_omega_n[i][1] - DH_t*(C_psi[i][1] + C_rho_x[i][1])/3;
			C_rho_n[i][0]   = C_rho_n[i][0]   - C_phi[i][0]*DH_t/3;
			C_rho_n[i][1]   = C_rho_n[i][1]   - C_phi[i][1]*DH_t/3;
		}

		for(i = 0; i < M*NUM; i++){												///4th, get K4
			C_omega[i][0]   = C_omega_o[i][0] - DH_t*(C_psi[i][0] + C_rho_x[i][0]);
			C_omega[i][1]   = C_omega_o[i][1] - DH_t*(C_psi[i][1] + C_rho_x[i][1]);
			C_rho[i][0]     = C_rho_o[i][0]   - C_phi[i][0]*DH_t;
			C_rho[i][1]     = C_rho_o[i][1]   - C_phi[i][1]*DH_t;	
		}
		Value_u_v(C_omega, C_u, C_v);
		Convolution(C_phi, C_psi);
		for(i = 0; i < M*NUM; i++){
			C_omega_n[i][0] = C_omega_n[i][0] - DH_t*(C_psi[i][0] + C_rho_x[i][0])/6;	
			C_omega_n[i][1] = C_omega_n[i][1] - DH_t*(C_psi[i][1] + C_rho_x[i][1])/6;
			C_rho_n[i][0]   = C_rho_n[i][0]   - C_phi[i][0]*DH_t/6;
			C_rho_n[i][1]   = C_rho_n[i][1]   - C_phi[i][1]*DH_t/6;
		}

		for(i = 0; i < M*NUM; i++){
			C_omega_o[i][0] = C_omega_n[i][0];	
			C_omega_o[i][1] = C_omega_n[i][1];	
			C_rho_o[i][0]   = C_rho_n[i][0];
			C_rho_o[i][1]   = C_rho_n[i][1];
		}
		Step = Step + 1;
	}while(Step < Step_Max);                 /// condition for loop end, t < 3.16
	
	fftw_free(C_psi);
	fftw_free(C_phi);
}
